R Notebook

Parameters

filtered_data <- read.csv("data/data-2023-09-11 (2).csv", header = TRUE)

Warning: cannot open file 'data/data-2023-09-11 (2).csv': No such file or directoryError in file(file, "rt") : cannot open the connection

data filtering convertion

mydata_genind

/// GENIND OBJECT /////////

 // 698 individuals; 6 loci; 48 alleles; size: 195 Kb

 // Basic content
   @tab:  698 x 48 matrix of allele counts
   @loc.n.all: number of alleles per locus (range: 5-12)
   @loc.fac: locus factor for the 48 columns of @tab
   @all.names: list of allele names for each locus
   @ploidy: ploidy of each individual  (range: 2-2)
   @type:  codom
   @call: df2genind(X = filtered_data[, 6:11], sep = "/", ncode = 6, ind.names = filtered_data$indv, 
    pop = filtered_data$Population, NA.char = "0/0", ploidy = 2, 
    type = "codom", strata = NULL, hierarchy = NULL)

 // Optional content
   @pop: population of each individual (group size range: 24-166)

Run basic.stats and render the result

library("pegas")
library(dplyr)
library(tibble)

result <- basic.stats(mydata_hierfstat)
df_resutl_basic<-as.data.frame(result$perloc)

# Weir and Cockrham estimates of Fstatistics - FIS and FST 
result_f_stats <- Fst(as.loci(mydata_genind))
result_f_stats <- result_f_stats[,2:3]
colnames(result_f_stats) <- c("Fst (W&C)", "Fis (W&C)")
result_f_stats <- merge(result_f_stats, df_resutl_basic, by="row.names",all.x=TRUE)
colnames(result_f_stats)[10] <- "Fst (Nei)"
colnames(result_f_stats)[12] <- "Fis (Nei)"
result_f_stats <- result_f_stats %>% column_to_rownames(., var = 'Row.names')
result_f_stats_selec <- result_f_stats %>% select(all_of(selected_stats))
result_f_stats_selec

G-statistic

library(hrbrthemes)


# compute the Gstats
result_f_stats <- result_f_stats %>% mutate(GST = 1-Hs/Ht)
result_f_stats <- result_f_stats %>% mutate("GST''" = (n_pop*(Ht-Hs))/((n_pop*Hs-Ht)*(1-Hs)))
result_f_stats

# Plot with linear trend
p2 <- ggplot(result_f_stats, aes(x=GST, y=Hs)) +
  geom_point() +
  geom_smooth(method=lm , color="red", se=FALSE) +
  theme_ipsum()
p2

Missing data

HW - Panmixia

library("pegas")

hw.test(as.loci(mydata_genind), B = 1000)

Linkage Disequilibrium

LD(as.loci(mydata_genind$tab), locus = c(1, 2), details = TRUE)

a<-LDscan(as.loci(mydata_hierfstat[,2:7]))

LDmap(a)


head(mydata_hierfstat[,2:7])


library("poppr")
pair.ia(mydata_genind, sample = 9)

shuffle df

#shuffled_matrices <- replicate(n_rep, mat[sample(nrow(mat)), ], simplify = FALSE)
shuffled_matrices <- replicate(n_rep, mat[sample(length(mat), replace = FALSE)], simplify = FALSE)
##################
# shuffle only the genotype and add the pop column later for each matrices. 
#in the loop? 
###############


# Create a list to store the wc
fst_df <- numeric(sequence_length)
fis_df <- numeric(sequence_length)

# Calculate the average for each shuffled matrix

# Iterate through the shuffled matrices
for (i in 1:n_rep) {
  # Calculate the statistics for the i-th matrix
  #HERE THE COLUMN POP
  merged_df <- cbind(level_pop, shuffled_matrices[[i]])
  result_f_stats <- wc(shuffled_matrices[[i]]) 
  result_f_stats <- as.data.frame(result_f_stats$per.loc)
  # Extract FST and FIS values
  fst_values <- result_f_stats$FST
  fis_values <- result_f_stats$FIS
  print( fst_values)
  # Assign values to the data frames
  fst_df <- cbind(fst_df, result_f_stats$FST)
  fis_df <- cbind(fis_df, result_f_stats$FIS)
}

# Set row names as in result_f_stats

rownames(fst_df) <- rownames(fis_df) <- rownames(result_f_stats)
result_FST <- fst_df[, -1]
fis_df <-fis_df[, -1]
vec <- seq(1, n_rep)
colnames(result_FST) <- colnames(fis_df) <- vec


result_FST[1,] 

count (result_f_stats[,1][1] > result_FST[1,] )
count <- sum(result_f_stats[,1][1] > result_FST[1, ])



# Initialize an empty data frame to store the counts
count_df <- data.frame(
  Greater = numeric(length(result_FST)),
  Smaller = numeric(length(result_FST))
)

# Compare the values in result_f_stats[1] to result_FST for each column
for (col in colnames(result_FST)) {
  greater_count <- sum(result_f_stats[1] > result_FST[, col])
  smaller_count <- sum(result_f_stats[1] < result_FST[, col])
  count_df$Greater[col] <- greater_count
  count_df$Smaller[col] <- smaller_count
}

# Print the count data frame
print(count_df)

######################## ########################